Context

column

The business situation

Our aluminum recyling company just bought a renewable energy manufacturing and services company to expand its earnings opportunities. The renewables devision dabbles in wind, clean energy technologies (very similar to the aluminum recycling clean technologies), and solar, a very new field for the company. The CFO would like to measure the impact of this new market on the earnings of the renewables division. To do this she commissions a project team.

Complications

  • The company has experienced very high volatility of earnings in aluminum and input energy markets. There is increased competition in these and their substitute and complementary markets. New technology is outpacing the efficiency and cost per unit of existing facilities in the company’s fleet of plants. The company continues to violate levels of risk tolerance and thresholds for total return set by the board of directors,

Key questions

Your CFO has a few questions for us:

  1. How do we characterize renewables variability and the impact of one market on another?

  2. What are the best combinations of renewables drivers?

  3. How much capital is needed to support a renewables earnings stream?

  4. How should the company plan to meet risk tolerances and thresholds for losses?

Data and workflow

For the renewables sector we select exchange traded funds (ETF) fromm the global renewables sector:

  • TAN for solar,

  • ICLN for clean technologies, and

  • PBW for wind.

These funds act as indices to effectively summarize the inputs, process, management, decisions, and outputs of various aspects of the renewables sector. Examining and analyzing this series will go a long way to helping the CFO understand the riskiness of these markets.

We load historical data on three ETFs, tranform prices into returns, and then further transform the returns into within-month correlations and standard deviations.

Our process includes

  • Review the stylized facts of volatility and relationships among three repesentative markets.

  • Develop market risk measures for each driver of earnings.

  • Apply corporate risk tolerances and thresholds to determine optimal collateral positions for each driver of earnings.

  • Determine optimal combinations of the drivers for maximum excess portfolio return relative to portfolio risk as well as the minimization of risk

  • Given the optimal maximum excess return per risk portfolio, determine the probable range of collateral needed to satisfy corporate risk tolerance and thresholds.

Returns

column

Renewables

TAN

ICLN

PBW

Return persistence

Volatility

column

TAN

PBW

ICLN

column

TAN-ICLN

Notes:

  • Interesting

  • Very interesting

TAN-PBW

ICLN-PBW

TAN-ICLN market spillover

TAN-PBW market spillover

ICLN-PBW market spillover

Loss

column

TAN

ICLN

PBW

Allocation

column

Efficient frontier

Sharpe mean CI

Sharpe standard deviation CI

Max Sharpe Ratio loss thresholds

Risky capital

Collateral

Summary

column

Minnie: TAN

  • TAN by it self

  • TAN-ICLN

  • TAN-PBW

Donald: ICLN

  • ICLN

Goofy: PBW

  • PBW

References

---
title: "Workbook: All in"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
    source_code: embed
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, eval = TRUE, warning=FALSE, message=FALSE)
library(flexdashboard)
#
options(digits = 4, scipen = 999999)
library(psych)
library(ggplot2)
library(GGally)
library(lubridate)
library(dplyr)
library(quantreg)
library(forecast)
library(tidyquant)
library(timetk)
library(quantmod)
library(matrixStats)
library(plotly)
library(quadprog)
#
#stocks_env <- new.env()
symbols <- c("TAN", "ICLN", "PBW") #c("ENE", "REP", "")

price_tbl <- tq_get(symbols) %>% select(date, symbol, price = adjusted)
# long format ("TIDY") price tibble for possible other work
return_tbl <- price_tbl %>% group_by(symbol) %>% tq_transmute(mutate_fun = periodReturn, period = "daily", type = "log", col_rename = "daily_return") %>% mutate(abs_return = abs(daily_return))
#str(return_tbl)
r_2 <- return_tbl %>% select(symbol, date, daily_return) %>% spread(symbol, daily_return)
r_2 <- xts(r_2, r_2$date)[-1, ]
storage.mode(r_2) <- "numeric"
r_2 <- r_2[, -1]
r_corr <- apply.monthly(r_2, FUN = cor)[,c(2, 3, 6)]
colnames(r_corr) <- c("TAN_ICLN", "TAN_PBW", "ICLN_PBW")
r_vols <- apply.monthly(r_2, FUN = colSds)
# 
corr_tbl <- r_corr %>% as_tibble() %>% mutate(date = index(r_corr)) %>% gather(key = assets, value = corr, -date)

vols_tbl <- r_vols %>% as_tibble() %>% mutate(date = index(r_vols)) %>% gather(key = assets, value = vols, -date) 
#
corr_vols <- merge(r_corr, r_vols)
corr_vols_tbl <- corr_vols %>% as_tibble() %>% mutate(date = index(corr_vols))
#
n <-  10000 # lots of trials, each a "day" or an "hour"
z <- rt(n, df = 30)
garch_sim_t <- function(n = 1000, df = 30, omega = 0.1, alpha = 0.8, phi = 0.05, mu = 0.01){
  n <- n # lots of trials, each a "day" or an "hour"
  # set.seed(seed)
  z <- rt(n, df = df) 
  e <-  z # store variates
  y <-  z # returns: store again in a different place
  sig2 <-  z^2 # create volatility series
  omega <-  omega #base variance
  alpha <-  alpha #vols Markov dependence on previous variance
  phi <-  phi # returns Markov dependence on previous period
  mu <-  mu # average return
  for (t in 2:n) { # Because of lag start at second
    e[t] <- sqrt(sig2[t])*z[t]           # 1. e is conditional on sig
    y[t] <-  mu + phi*(y[t-1]-mu) + e[t] # 2. generate returns
    sig2[t+1] <-  omega + alpha * e[t]^2 # 3. generate new sigma^2
    }
  return <- list(
    sim_df_vbl <- data_frame(t = 1:n, z = z, y = y, e = e, sig = sqrt(sig2)[-(n+1)] ),
    sim_df_title <- data_frame(t = 1:n, "1. Unconditional innovations" = z, "4. Conditional returns" = y, "3. Conditional innovations" = e, "2. Conditional volatility" = sqrt(sig2)[-(n+1)] )
  )
}
#
# convert prices from tibble to xts
price_etf <- price_tbl %>% spread(symbol, price)
price_etf <- xts(price_etf, price_etf$date)
storage.mode(price_etf) <- "numeric" #select(TAN, ICLN, PBW) # 3 risk factors (rf)
price_etf <- price_etf[, -1]
price_0 <- as.numeric(tail(price_etf, 1))
shares <- c(60000, 75000, 50000)
#price_last <- price_etf[length(price_etf$TAN), 3:5] #(TAN, ICLN, PBW) %>% as.vector()
w <- as.numeric(shares * price_0)
return_hist <- r_2
# Fan these across the length and breadth of the risk factor series
weights_rf <- matrix(w, nrow=nrow(return_hist), ncol=ncol(return_hist), byrow=TRUE)
## We need to compute exp(x) - 1 for very small x: expm1 accomplishes this
loss_rf <- -rowSums(expm1(return_hist) * weights_rf)
loss_df <- data_frame(loss = loss_rf, distribution = rep("historical", each = length(loss_rf)))
#
ES_calc <- function(data, prob){
  threshold <- quantile(data, prob)
  result <- mean(data[data > threshold])
}
#
n_sim <- 1000
n_sample <- 100
prob <- 0.95
ES_sim <- replicate(n_sim, ES_calc(sample(loss_rf, n_sample, replace = TRUE), prob))
summary(ES_sim)
#
#summary(ES_sim)
#
# mean excess plot to determine thresholds for extreme event management
data <- as.vector(loss_rf) # data is purely numeric
umin <-  min(data)         # threshold u min
umax <-  max(data) - 0.1   # threshold u max
nint <- 100                # grid length to generate mean excess plot
grid_0 <- numeric(nint)    # grid store
e <- grid_0                # store mean exceedances e
upper <- grid_0            # store upper confidence interval
lower <- grid_0            # store lower confidence interval
u <- seq(umin, umax, length = nint) # threshold u grid
alpha <- 0.95                  # confidence level
for (i in 1:nint) {
    data <- data[data > u[i]]  # subset data above thresholds
    e[i] <- mean(data - u[i])  # calculate mean excess of threshold
    sdev <- sqrt(var(data))    # standard deviation
    n <- length(data)          # sample size of subsetted data above thresholds
    upper[i] <- e[i] + (qnorm((1 + alpha)/2) * sdev)/sqrt(n) # upper confidence interval
    lower[i] <- e[i] - (qnorm((1 + alpha)/2) * sdev)/sqrt(n) # lower confidence interval
  }
mep_df <- data.frame(threshold = u, threshold_exceedances = e, lower = lower, upper = upper)
loss_excess <- loss_rf[loss_rf > u] - u
quantInv <- function(distr, value) ecdf(distr)(value)
u_prob <- quantInv(loss_rf, 200000)
ES_mep <- mean(loss_rf[loss_rf > u_prob])
stat_fun <- function(x, na.rm = TRUE, ...) {
  library(moments)
    # x     = numeric vector
    # na.rm = boolean, whether or not to remove NA's
    # ...   = additional args passed to quantile
    c(mean     = mean(x, na.rm = na.rm),
      stdev    = sd(x, na.rm = na.rm),
      skewness = skewness(x, na.rm = na.rm),
      kurtosis = kurtosis(x, na.rm = na.rm),
      quantile(x, na.rm = na.rm, ...)) 
}
#
contract <- 1 # billion
working <- 0.100 # billion
sigma_wc <- 0.025 # billion
sigma <- 0.25
threshold <- -0.12 # percentage return
alpha <- 0.05 # tolerance
risky <- 0.1 # percentage return on the risky asset
riskless <- 0.02 # time value of cash -- no risk
z_star <- qnorm(alpha)
w <- (threshold-riskless) / (risky - riskless + sigma*z_star)# Tukey-Box-Hunter fence analysis of outliers
#
k <- 1:20 # days in a business month
col_names <- paste0("lag_", k)
#
# remove abs_return the fourth column
return_lags <- return_tbl[, -4] %>%
  tq_mutate(
  select     = daily_return,
  mutate_fun = lag.xts,
  k          = k,
  col_rename = col_names
  )
return_autocors <- return_lags %>%
  gather(key = "lag", value = "lag_value", -c(symbol, date, daily_return)) %>%
  mutate(lag = str_sub(lag, start = 5) %>% as.numeric) %>%
  group_by(symbol, lag) %>%
  summarize(
    cor = cor(x = daily_return, y = lag_value, use = "pairwise.complete.obs"),
    upper_95 = 2/(n())^0.5,
    lower_95 = -2/(n())^0.5
  )
return_absautocors <- return_autocors %>%
  ungroup() %>%
  mutate(
    lag = as_factor(as.character(lag)),
    cor_abs = abs(cor)
  ) %>%
  select(lag, cor_abs) %>%
  group_by(lag)
#
## INPUTS: r vector
## OUTPUTS: list of scalars (mean, sd, median, skewness, kurtosis)
data_moments <- function(data){
  library(moments)
  library(matrixStats)
  mean <- colMeans(data)
  median <- colMedians(data)
  sd <- colSds(data)
  IQR <- colIQRs(data)
  skewness <- skewness(data)
  kurtosis <- kurtosis(data)
  result <- data.frame(mean = mean, median = median, std_dev = sd, IQR = IQR, skewness = skewness, kurtosis = kurtosis)
  return(result)
}
#
#
port_sample <- function(return, n_sample = 252, stat = "mean")
{
  R <-  return # daily returns
  n <- dim(R)[1]
  N <- dim(R)[2]
  R_boot <-  R[sample(1:n, n_sample),] # sample returns
  r_free <- 0.03 / 252 # daily
  mean_vect <-  apply(R_boot,2,mean)
  cov_mat <-  cov(R_boot)
  sd_vect <-  sqrt(diag(cov_mat))
  A_mat <-  cbind(rep(1,N),mean_vect) 
  mu_P <-  seq(-.01,.01,length=300)                              
  sigma_P <-  mu_P 
  weights <-  matrix(0,nrow=300,ncol=N) 
  for (i in 1:length(mu_P))  
  {
    b_vec <-  c(1,mu_P[i])  
    result <-  
      solve.QP(Dmat=2*cov_mat,dvec=rep(0,N),Amat=A_mat,bvec=b_vec,meq=2)
    sigma_P[i] <-  sqrt(result$value)
    weights[i,] <-  result$solution
  }
  sharpe <- (mu_P - r_free)/sigma_P ## compute Sharpe's ratios
  ind_max <-  (sharpe == max(sharpe)) ## Find maximum Sharpe's ratio
  ind_min <-  (sigma_P == min(sigma_P)) ## find the minimum variance portfolio
  ind_eff <-  (mu_P > mu_P[ind_min]) ## finally the efficient fr(aes(x = 0, y = r_free), colour = "red")ontier
  result <- switch(stat,
    "mean"  = mu_P[ind_max],
    "sd"    = sigma_P[ind_max]
    )
  return(result)
}
#
```

Context
==============================================

column{.tabset}
------------------------------------------------------

### The business situation

Our aluminum recyling company just bought a renewable energy manufacturing and services company to expand its earnings opportunities. The renewables devision dabbles in wind, clean energy technologies (very similar to the aluminum recycling clean technologies), and solar, a very new field for the company. The CFO would like to measure the  impact of this new market on the earnings of the renewables division. To do this she commissions a project team.

**Complications**

- The company has experienced very high volatility of earnings in aluminum and input energy markets. There is increased competition in these and their substitute and complementary markets. New technology is outpacing the efficiency and cost per unit of existing facilities in the company's fleet of plants. The company continues to violate levels of risk tolerance and thresholds for total return set by the board of directors,


***Key questions***

Your CFO has a few questions for us:

1. How do we characterize renewables variability and the impact of one market on another? 

2. What are the best combinations of renewables drivers?

3. How much capital is needed to support a renewables earnings stream?

4. How should the company plan to meet risk tolerances and thresholds for losses?

### Data and workflow

For the renewables sector we select [exchange traded funds (ETF)](https://www.investopedia.com/terms/e/etf.asp) fromm the [global renewables sector](https://www.etf.com/channels/renewable-energy-etfs): 

- TAN for solar, 

- ICLN for clean technologies, and 

- PBW for wind. 

These funds act as indices to effectively summarize the inputs, process, management, decisions, and outputs of various aspects of the renewables sector. Examining and analyzing this series will go a long way to helping the CFO understand the riskiness of these markets. 

We load historical data on three ETFs, tranform prices into returns, and then further transform the returns into within-month correlations and standard deviations.

Our process includes 

- Review the stylized facts of volatility and relationships among three repesentative markets.

- Develop market risk measures for each driver of earnings.

- Apply corporate risk tolerances and thresholds to determine optimal collateral positions for each driver of earnings.

- Determine optimal combinations of the drivers for maximum excess portfolio return relative to portfolio risk as well as the minimization of risk

- Given the optimal maximum excess return per risk portfolio, determine the probable range of collateral needed to satisfy corporate risk tolerance and thresholds.

Returns
=======================================================================

column {.sidebar}
----------------------------------------------------

Summary

#### Overall Notes #### Persistence Notes #### Outliers Notes column {.tabset} ---------------------------------------------------- ### Renewables ```{r moments} return_plot <- return_tbl %>% select(date, symbol, daily_return) %>% spread(symbol, daily_return) ggpairs(return_plot) ``` ### TAN ```{r} ggtsdisplay(return_plot$TAN, plot.type = "histogram", main = "TAN daily returns") ``` ### ICLN ```{r} ggtsdisplay(return_plot$ICLN, plot.type = "histogram", main = "ICLN daily returns") ``` ### PBW ```{r} ggtsdisplay(return_plot$PBW, plot.type = "histogram", main = "PBW daily returns") ``` ### Return persistence ```{r plotreturnabscors, exercise = TRUE} # Tukey's fence upper_bound <- 1.5*IQR(return_absautocors$cor_abs) %>% signif(3) p <- return_absautocors %>% ggplot(aes(x = fct_reorder(lag, cor_abs, .desc = TRUE) , y = cor_abs)) + # Add boxplot geom_boxplot(color = palette_light()[[1]]) + # Add horizontal line at outlier break point geom_hline(yintercept = upper_bound, color = "red") + annotate("text", label = paste0("Outlier threshold = ", upper_bound), x = 24.5, y = upper_bound + .03, color = "red") + # Aesthetics expand_limits(y = c(0, 0.5)) + theme_tq() + labs( title = paste0("Absolute Autocorrelations: Lags ", rlang::expr_text(k)), x = "Lags" ) + theme( legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1) ) ggplotly(p) ``` Volatility ==================================================== column {.sidebar} ----------------------------------------------------

Summary

### Monthly volatility Notes ### Monthly Correlation Notes column {.tabset} ---------------------------------------------------- ### TAN ```{r} ggtsdisplay(r_vols$TAN, plot.type = "histogram", main = "TAN monthly volatility") ``` ### PBW ```{r} ggtsdisplay(r_vols$ICLN, plot.type = "histogram", main = "ICLN monthly volatility") ``` ### ICLN ```{r} ggtsdisplay(r_vols$PBW, plot.type = "histogram", main = "PBW monthly volatility") ``` column {.tabset} ---------------------------------------------------- ### TAN-ICLN ```{r} ggtsdisplay(r_corr$TAN_ICLN, plot.type = "histogram", main = "TAN-ICLN monthly correlation") ``` Notes: - Interesting - Very interesting ### TAN-PBW ```{r} ggtsdisplay(r_corr$TAN_PBW, plot.type = "histogram", main = "TAN-PBW monthly correlation") ``` ### ICLN-PBW ```{r} ggtsdisplay(r_corr$ICLN_PBW, plot.type = "histogram", main = "ICLN-PBW monthly correlation") ``` ### TAN-ICLN market spillover ```{r rqplot-TAN-ICLN} p <- ggplot(corr_vols_tbl, aes(x = ICLN, y = TAN_ICLN)) + geom_point() + ggtitle("TAN-ICLN Interaction") + geom_quantile(quantiles = c(0.10, 0.90)) + geom_quantile(quantiles = 0.5, linetype = "longdash") + geom_density_2d(colour = "red") ggplotly(p) ``` ### TAN-PBW market spillover ```{r rqplot-TAN-PBW} p <- ggplot(corr_vols_tbl, aes(x = PBW, y = TAN_PBW)) + geom_point() + ggtitle("TAN-PBW Interaction") + geom_quantile(quantiles = c(0.10, 0.90)) + geom_quantile(quantiles = 0.5, linetype = "longdash") + geom_density_2d(colour = "red") ggplotly(p) ``` ### ICLN-PBW market spillover ```{r rqplot-ICLN-PBW} p <- ggplot(corr_vols_tbl, aes(x = PBW, y = ICLN_PBW)) + geom_point() + ggtitle("ICLN-PBW Interaction") + geom_quantile(quantiles = c(0.10, 0.90)) + geom_quantile(quantiles = 0.5, linetype = "longdash") + geom_density_2d(colour = "red") ggplotly(p) ``` Loss ============================================ column {.sidebar} ----------------------------------------------------

Pure plays

### TAN only Notes ### ICLN only Notes ### PBW only Notes column {.tabset} -------------------------------------------- ### TAN ```{r TANloss} # shares <- c(-215000, 0, 0) price_last <- c(1, 0, 0) * price_0 #(TAN, ICLN, PBW) %>% as.vector() w <- as.numeric(shares * price_last) return_hist <- r_2 # Fan these across the length and breadth of the risk factor series weights_rf <- matrix(w, nrow=nrow(return_hist), ncol=ncol(return_hist), byrow=TRUE) ## We need to compute exp(x) - 1 for very small x: expm1 accomplishes this loss_rf <- -rowSums(expm1(return_hist) * weights_rf) loss_df <- data_frame(loss = loss_rf, distribution = rep("historical", each = length(loss_rf))) # ES_calc <- function(data, prob){ threshold <- quantile(data, prob) result <- mean(data[data > threshold]) } # n_sim <- 1000 n_sample <- 100 prob <- 0.95 ES_sim <- replicate(n_sim, ES_calc(sample(loss_rf, n_sample, replace = TRUE), prob)) # sim <- ES_sim low <- quantile(sim, 0.025) high <- quantile(sim, 0.975) sim_df <- data_frame(sim = sim) title <- "TAN: Expected Shortfall simulation" p <- ggplot(data = sim_df, aes(x = sim)) p <- p + geom_histogram(binwidth = 1000, aes(y = 1000*(..density..)), alpha = 0.4) p <- p + ggtitle(title) p <- p + geom_vline(xintercept = low, color = "red", size = 1.5 ) + geom_vline(xintercept = high, color = "red", size = 1.5) p <- p + annotate("text", x = low, y = 0.01, label = paste("L = ", round(low, 2))) + annotate("text", x = high, y = 0.01, label = paste("U = ", round(high, 2))) + ylab("density") + xlab("expected shortfall") + theme_bw() ggplotly(p) ``` ### ICLN ```{r ICLNloss} # shares <- c(0, 284000, 0) price_last <- c(0, 1, 0) * price_0 w <- as.numeric(shares * price_last) return_hist <- r_2 # Fan these across the length and breadth of the risk factor series weights_rf <- matrix(w, nrow=nrow(return_hist), ncol=ncol(return_hist), byrow=TRUE) ## We need to compute exp(x) - 1 for very small x: expm1 accomplishes this loss_rf <- -rowSums(expm1(return_hist) * weights_rf) loss_df <- data_frame(loss = loss_rf, distribution = rep("historical", each = length(loss_rf))) # ES_calc <- function(data, prob){ threshold <- quantile(data, prob) result <- mean(data[data > threshold]) } # n_sim <- 1000 n_sample <- 100 prob <- 0.95 ES_sim <- replicate(n_sim, ES_calc(sample(loss_rf, n_sample, replace = TRUE), prob)) # sim <- ES_sim low <- quantile(sim, 0.025) high <- quantile(sim, 0.975) sim_df <- data_frame(sim = sim) title <- "ICLN: Expected Shortfall simulation" p <- ggplot(data = sim_df, aes(x = sim)) p <- p + geom_histogram(binwidth = 1000, aes(y = 1000*(..density..)), alpha = 0.4) p <- p + ggtitle(title) p <- p + geom_vline(xintercept = low, color = "red", size = 1.5 ) + geom_vline(xintercept = high, color = "red", size = 1.5) p <- p + annotate("text", x = low, y = 0.01, label = paste("L = ", round(low, 2))) + annotate("text", x = high, y = 0.01, label = paste("U = ", round(high, 2))) + ylab("density") + xlab("expected shortfall") + theme_bw() ggplotly(p) ``` ### PBW ```{r PBWloss} # shares <- c(0, 0, 12500) price_last <- c(0, 0, 1) * price_0 return_hist <- r_2 # Fan these across the length and breadth of the risk factor series weights_rf <- matrix(w, nrow=nrow(return_hist), ncol=ncol(return_hist), byrow=TRUE) ## We need to compute exp(x) - 1 for very small x: expm1 accomplishes this loss_rf <- -rowSums(expm1(return_hist) * weights_rf) loss_df <- data_frame(loss = loss_rf, distribution = rep("historical", each = length(loss_rf))) # ES_calc <- function(data, prob){ threshold <- quantile(data, prob) result <- mean(data[data > threshold]) } # n_sim <- 1000 n_sample <- 100 prob <- 0.95 ES_sim <- replicate(n_sim, ES_calc(sample(loss_rf, n_sample, replace = TRUE), prob)) # sim <- ES_sim low <- quantile(sim, 0.025) high <- quantile(sim, 0.975) sim_df <- data_frame(sim = sim) title <- "PBW: Expected Shortfall simulation" p <- ggplot(data = sim_df, aes(x = sim)) p <- p + geom_histogram(binwidth = 1000, aes(y = 1000*(..density..)), alpha = 0.4) p <- p + ggtitle(title) p <- p + geom_vline(xintercept = low, color = "red", size = 1.5 ) + geom_vline(xintercept = high, color = "red", size = 1.5) p <- p + annotate("text", x = low, y = 0.01, label = paste("L = ", round(low, 2))) + annotate("text", x = high, y = 0.01, label = paste("U = ", round(high, 2))) + ylab("density") + xlab("expected shortfall") + theme_bw() ggplotly(p) ``` Allocation ============================================ column {.sidebar} -------------------------------------------- ```{r eff-frontier-calc} R <- r_2 # daily returns n <- dim(R)[1] N <- dim(R)[2] R_boot <- R[sample(1:n, 252),] # sample returns r_free <- 0.03 / 252 # daily mean_vect <- apply(R_boot,2,mean) cov_mat <- cov(R_boot) sd_vect <- sqrt(diag(cov_mat)) A_mat <- cbind(rep(1,N),mean_vect) mu_P <- seq(-.01,.01,length=300) sigma_P <- mu_P weights <- matrix(0,nrow=300,ncol=N) for (i in 1:length(mu_P)) { b_vec <- c(1,mu_P[i]) result <- solve.QP(Dmat=2*cov_mat,dvec=rep(0,N),Amat=A_mat,bvec=b_vec,meq=2) sigma_P[i] <- sqrt(result$value) weights[i,] <- result$solution } # make a data frame of the mean and standard deviation results sigma_mu_df <- data_frame(sigma_P = sigma_P, mu_P = mu_P) names_R <- c("TAN", "ICLN", "PBW") # sharpe ratio and minimum variance portfolio analysis sharpe <- (mu_P - r_free)/sigma_P ## compute Sharpe's ratios ind_max <- (sharpe == max(sharpe)) ## Find maximum Sharpe's ratio ind_min <- (sigma_P == min(sigma_P)) ## find the minimum variance portfolio ind_eff <- (mu_P > mu_P[ind_min]) ## finally the efficient fr(aes(x = 0, y = r_free), colour = "red")ontier w_max <- weights[ind_max,] w_min <- weights[ind_min,] value <- 1000000 ``` ### Optimal weights We use these maximum Sharpe Ratio weights to form a risky asset loss series: TAN: `r round(w_max[1] * 100, 2)`%, ICLN: `r round(w_max[2] * 100, 2)`%, PBW: `r round(w_max[3] * 100, 2)`% of total portfolio value of USD `r value` with a daily return of `r round(mu_P[ind_max]*100, 2)`% and standard deviation of `r round(sigma_P[ind_max]*100, 2)`% These weights will produce a minimum variance portfolio: TAN: `r round(w_min[1] * 100, 2)`%, ICLN: `r round(w_min[2] * 100, 2)`%, PBW: `r round(w_min[3] * 100, 2)`% of total portfolio value of USD `r value` with a daily portfolio return of `r round(mu_P[ind_min]*100, 2)`% and standard deviation of `r round(sigma_P[ind_min]*100, 2)`% ### Thresholds - overall ### Loss range ### Collateral concerns How much collateral do we need to meet company risk tolerances and thresholds? column {.tabset} -------------------------------------------- ### Efficient frontier ```{r eff-frontier} col_P <- ifelse(mu_P > mu_P[ind_min], "blue", "grey") # discriminate efficient and inefficient portfolios sigma_mu_df$col_P <- col_P p <- ggplot(sigma_mu_df, aes(x = sigma_P, y = mu_P, group = 1)) p <- p + geom_line(aes(colour=col_P, group = col_P), size = 1.1) + scale_colour_identity() p <- p + geom_abline(intercept = r_free, slope = (mu_P[ind_max]-r_free)/sigma_P[ind_max], color = "red", size = 1.1) p <- p + geom_point(aes(x = sigma_P[ind_max], y = mu_P[ind_max]), color = "green", size = 4) p <- p + geom_point(aes(x = sigma_P[ind_min], y = mu_P[ind_min]), color = "red", size = 4) ## show min var portfolio #p ggplotly(p) ``` ### Sharpe mean CI ```{r sampledmean-ex} port_mean <- replicate(1000, port_sample(R, n_sample = 252, stat = "mean")) sim <- port_mean * 252 low <- quantile(sim, 0.025) high <- quantile(sim, 0.975) sim_df <- data_frame(sim = sim) title <- "Tangency portfolio sampled mean simulation" p <- ggplot(data = sim_df, aes(x = sim)) p <- p + geom_histogram(alpha = 0.7) p <- p + ggtitle(title) p <- p + geom_vline(xintercept = low, color = "red", size = 1.5 ) + geom_vline(xintercept = high, color = "red", size = 1.5) p <- p + annotate("text", x = low + 0.01, y = 200, label = paste("L = ", round(low, 2))) + annotate("text", x = high, y = 200, label = paste("U = ", round(high, 2))) + ylab("density") + xlab("daily mean: max Sharpe Ratio") + theme_bw() ggplotly(p) ``` ### Sharpe standard deviation CI ```{r sampledsd-ex} port_mean <- replicate(1000, port_sample(R, n_sample = 252, stat = "sd")) sim <- port_mean * 252 low <- quantile(sim, 0.025) high <- quantile(sim, 0.975) sim_df <- data_frame(sim = sim) title <- "Tangency portfolio sampled standard deviation simulation" p <- ggplot(data = sim_df, aes(x = sim)) p <- p + geom_histogram(alpha = 0.7) p <- p + ggtitle(title) p <- p + geom_vline(xintercept = low, color = "red", size = 1.5 ) + geom_vline(xintercept = high, color = "red", size = 1.5) p <- p + annotate("text", x = low + 0.1, y = 200, label = paste("L = ", round(low, 2))) + annotate("text", x = high, y = 200, label = paste("U = ", round(high, 2))) + ylab("density") + xlab("daily mean: max Sharpe Ratio") + theme_bw() ggplotly(p) ``` ### Max Sharpe Ratio loss thresholds ```{r mepcalc} price_last <- price_0 value <- 1000000 # portfolio value w_0 <- w_max # wwights -- e.g., min variance or max sharpe shares <- value * (w_0/price_last) w <- as.numeric(shares * price_last) return_hist <- r_2 # Fan these across the length and breadth of the risk factor series weights_rf <- matrix(w, nrow=nrow(return_hist), ncol=ncol(return_hist), byrow=TRUE) ## We need to compute exp(x) - 1 for very small x: expm1 accomplishes this loss_rf <- -rowSums(expm1(return_hist) * weights_rf) loss_df <- data_frame(loss = loss_rf, distribution = rep("historical", each = length(loss_rf))) # data <- as.vector(loss_rf[loss_rf > 0]) # data is purely numeric umin <- min(data) # threshold u min umax <- max(data) - 0.1 # threshold u max nint <- 100 # grid length to generate mean excess plot grid_0 <- numeric(nint) # grid store e <- grid_0 # store mean exceedances e upper <- grid_0 # store upper confidence interval lower <- grid_0 # store lower confidence interval u <- seq(umin, umax, length = nint) # threshold u grid alpha <- 0.95 # confidence level for (i in 1:nint) { data <- data[data > u[i]] # subset data above thresholds e[i] <- mean(data - u[i]) # calculate mean excess of threshold sdev <- sqrt(var(data)) # standard deviation n <- length(data) # sample size of subsetted data above thresholds upper[i] <- e[i] + (qnorm((1 + alpha)/2) * sdev)/sqrt(n) # upper confidence interval lower[i] <- e[i] - (qnorm((1 + alpha)/2) * sdev)/sqrt(n) # lower confidence interval } mep_df <- data.frame(threshold = u, threshold_exceedances = e, lower = lower, upper = upper) ``` ```{r loss-mep} # Voila the plot => you may need to tweak these limits! p <- ggplot(mep_df, aes( x= threshold, y = threshold_exceedances)) + geom_line() + geom_line(aes(x = threshold, y = lower), colour = "red") + geom_line(aes(x = threshold, y = upper), colour = "red") + annotate("text", x = mean(mep_df$threshold), y = max(mep_df$upper)+100, label = "upper 95%") + annotate("text", x = mean(mep_df$threshold), y = min(mep_df$lower) - 100, label = "lower 5%") + ggtitle("Mean Excess Plot: maximum Sharpe Ratio portfolio") + ylab("threshold exceedances") ggplotly(p) ``` ### Risky capital ```{r loss-capital} n_sim <- 1000 n_sample <- 100 prob <- 0.95 ES_sim <- replicate(n_sim, ES_calc(sample(loss_rf, n_sample, replace = TRUE), prob)) # sim <- ES_sim low <- quantile(sim, 0.025) high <- quantile(sim, 0.975) sim_df <- data_frame(sim = sim) title <- paste0("Loss Capital Simulation: alpha = ", alpha*100, "% bounds") p <- ggplot(data = sim_df, aes(x = sim)) p <- p + geom_histogram(alpha = 0.4) p <- p + ggtitle(title) p <- p + geom_vline(xintercept = low, color = "red", size = 1.5 ) + geom_vline(xintercept = high, color = "red", size = 1.5) p <- p + annotate("text", x = low, y = 100, label = paste("L = ", round(low, 2))) + annotate("text", x = high, y = 100, label = paste("U = ", round(high, 2))) + ylab("density") + xlab("expected shortfall") + theme_bw() ggplotly(p) ``` ### Collateral ```{r collateral} options(digits = 2, scipen = 99999) # r_f <- 0.03 # per annum mu <- mu_P[ind_max] * 252 # pull mu and sigma for tangency portfolio sigma <- sigma_P[ind_max] * sqrt(252) # sigma_p <- seq(0, sigma + 0.25, length.out = 100) mu_p <- r_f + (mu - r_f)*sigma_p/sigma w <- sigma_p / sigma threshold <- -0.12 alpha <- 0.05 z_star <- qnorm(alpha) w_star <- (threshold-r_f) / (mu - r_f + sigma*z_star) sim_df <- data_frame(sigma_p = sigma_p, mu_p = mu_p, w = w) # label_42 <- paste(alpha*100, "% alpha, ", threshold*100, "% threshold, \n", round(w_star*100, 2), "% risky asset", sep = "") label_0 <- paste(0*100, "% risky asset", sep = "") label_100 <- paste(1.00*100, "% risky asset", sep = "") p <- ggplot(sim_df, aes(x = sigma_p, y = mu_p)) + geom_line(color = "blue", size = 1.1) + geom_point(aes(x = 0.0 * sigma, y = r_f + (mu-r_f)*0.0), color = "red", size = 3.0) + annotate("text", x = 0.2 * sigma, y = r_f + (mu-r_f)*0.0 + 0.01, label = label_0) + geom_point(aes(x = w_star * sigma, y = r_f + (mu-r_f)*w_star), shape = 21, color = "red", fill = "white", size = 4, stroke = 4) + annotate("text", x = w_star * sigma + .2, y = r_f + (mu-r_f)*w_star + 0.1, label = label_42) + geom_point(aes(x = 1.0 * sigma, y = r_f + (mu-r_f)*1.00), color = "red", size = 3.0) + annotate("text", x = 1.0 * sigma, y = r_f + (mu-r_f)*1.00 + 0.01, label = label_100) + xlab("standard deviation of portfolio return") + ylab("mean of portfolio return") + ggtitle("Risk-return tradeoff of cash and risky asset") ggplotly(p) ``` Summary ============================================ column {.sidebar} -------------------------------------------- ### Overall Notes column {.tabset} -------------------------------------------- ### Minnie: TAN - TAN by it self - TAN-ICLN - TAN-PBW ### Donald: ICLN - ICLN ### Goofy: PBW - PBW References ============================================